# Dickstein, Ho, and Mark (2023)
# This file defines the cost variables needed in the premium setting equation.

# Preliminaries
setwd("../../library")
source("PreliminariesCode.R")
library(readr)
library(knitr)
library(readxl)
library(numDeriv)
library(yaml)
library(haven)
try(library(Hmisc))
library(EnvStats)

run <- "lowrisk_simplemh_censor0.2_main"

main <- function() {
  # Read data
  data_path <- paste0(project_folder, "/analysis/estimationcode/specs/", run, "/output/diag_exploded.rda")
  load(data_path)
  
  # Predict the costs
  predict_cost(4197169, param_data)
}

predict_cost <- function(seed, param_data) {
  set.seed(seed)
  param_data <- param_data %>% 
    mutate(
      lambda = rexp(n(), rate=alpha_i),
      simcost=(x_av+omega_i*x_av^2)*lambda,
      wsimcost=simcost*(simcost < 89.37073) + 89.37073*(simcost > 89.37073),
      cbar = (x_av + omega_i * x_av^2)/alpha_i,
      w_cbar = cbar*(cbar < 89.37073) + 89.37073*(cbar > 89.37073),
      w95_cbar = cbar*(cbar < 47.32818) + 47.32818*(cbar > 47.32818),
      w80_cbar = cbar*(cbar < 20) + 20*(cbar > 20),
      nominal_av = 0.6*(metal=="2")+0.7*(metal=="3")+0.8*(metal=="4")+0*(1-insurance),
      kappa_ratio = nominal_av/x_av,
      weighted_simcost = kappa_ratio * simcost, 
      weighted_wsimcost = kappa_ratio * wsimcost, 
      weighted_cbar = kappa_ratio * cbar,
      weighted_wcbar = kappa_ratio * w_cbar,
      weighted95_cbar = kappa_ratio * w95_cbar,
      weighted80_cbar = kappa_ratio * w80_cbar,
      grossp = grossprem/100,
      s_simcost = s_ij * weighted_simcost,
      s_wsimcost = s_ij * weighted_wsimcost,
      s_cbar = s_ij * weighted_cbar,
      s_wcbar = s_ij * weighted_wcbar,
      s_95cbar = s_ij * weighted95_cbar,
      s_80cbar = s_ij * weighted80_cbar,
      s_p = s_ij * p,
      s_cost = s_ij * cost,
      s_grossp = s_ij * grossp)   

  grouped <- param_data %>% 
    group_by(
      constructed_plan_year, PAYER_ID, insurance, best_guess_ra, year,
      nominal_av, x_av, kappa_ratio) %>%
      summarise(
        mean_s_ij = mean(s_ij),
        sum_s_ij = sum(s_ij),
        p=mean(p),
        cost=mean(cost),
        grossp=mean(grossp),
        cbar=mean(cbar),
        w_cbar=mean(w_cbar),
        w95_cbar=mean(w95_cbar),
        w80_cbar=mean(w80_cbar),
        s_simcost=sum(s_simcost),
        s_wsimcost=sum(s_wsimcost),
        s_cbar=sum(s_cbar),
        s_wcbar=sum(s_wcbar),
        s_95cbar=sum(s_95cbar),
        s_80cbar=sum(s_80cbar),
        s_p=sum(s_p),
        s_cost=sum(s_cost),
        s_grossp=sum(s_grossp),
        weighted_cbar=mean(weighted_cbar),
        weighted_wcbar=mean(weighted_wcbar),
        subscriberid=sum(choice),
        num_subsidized=sum(subsidized*choice))

  output_path <- paste0(project_folder, "/analysis/prem_setting/output")
  dir.create(output_path, recursive = TRUE)
  write.dta(grouped, file=paste0(output_path, "/predicted_prem_", seed, ".dta"))
}

main()
